function F = solve_eqm_cf1_inner(param, const, integ, eqm)
%% Counterfactual 1: Baseline

    % fix P_s = 1, e=1.
    P_s     = param.P_s;
    e       = param.e;     

    % read in parameters
%     kappa   = param.kappa;
    alpha_m = param.alpha_m; 
    alpha_s = param.alpha_s;
%     beta_v  = param.beta_v; 
%     beta_m  = param.beta_m;    
    sigma   = param.sigma; 
%     f_vH    = param.f_vH;
%     f_vF    = param.f_vF;
    mu_E    = param.mu_E;
    sigma_E = param.sigma_E;   
    w       = param.w;
    Q       = param.Q;
    Nf      = param.Nf;
    
    % read in constants    
    gamma   = const.gamma;
    k0      = const.k0;
%     C_m     = const.C_m;
%     C_v0    = const.C_v0;
%     C_x0    = const.C_x0;   
%     phi_v   = const.phi_v;
    phi_y   = const.phi_y;
    phi_s   = const.phi_s;
    D_F     = const.D_F;
    lnzQ    = const.lnzQ;
        
    % read in all the integral from previous round
%     E_zm0 = integ.E_zm0;
%     E_zm1 = integ.E_zm1;
%     E_zv0 = integ.E_zv0;
%     E_zv1 = integ.E_zv1;
    E_zx0 = integ.E_zx0;
    E_zx1 = integ.E_zx1;
        
    % adjust for level in case integrals too large
    rescale = max(1,10^floor(log10(max(E_zx1))-1));
    E_zx0 = E_zx0/rescale;
    E_zx1 = E_zx1/rescale;
%     E_zv0 = E_zv0/(rescale^(1/beta_v));
%     E_zv1 = E_zv1/(rescale^(1/beta_v));
%     E_zm0 = E_zm0/(rescale^(1/beta_m));        
%     E_zm1 = E_zm1/(rescale^(1/beta_m));      

    
    newmbar = eqm.newmbar;


%% Solve the equilibrium    

    % Setting
    tol = 1e-12;
    maxIter = 1000;

    % Initial guess 
    D_Hnew = eqm.D_H;
    c_new = eqm.c;    
    
    % Iteration
    err = 1;
    iter = 0;
    while err >= tol && iter < maxIter

        % new guess
        D_H = D_Hnew;
        c = c_new;

        % D(q,E)
        D0 = D_H;
        D1 = D_H + e^sigma.*D_F;       
        
        % Pi(q)
        C_x = (sigma/(sigma-1)*w.^(1-alpha_m-alpha_s)).^(1-sigma);
        Pi0 = D0.*(c.^alpha_m*P_s^alpha_s).^(1-sigma).*C_x;
        Pi1 = D1.*(c.^alpha_m*P_s^alpha_s).^(1-sigma).*C_x;     
               
        % X(q,E)
        X0 = Pi0.*E_zx0;
        X1 = Pi1.*E_zx1;
        X_m = alpha_m*(sigma-1)/sigma*(X0+X1);

        % P(q)^(1-sigma)
        P0sig = X0./D0;
        P1sig = X1./D1;
        Psig = X0./D0 + X1./D1;
        
        % D_m(q)
        TempD = c.^(sigma-1).*X_m;
        TempD(isnan(TempD))=0;
        D_m = newmbar*sum(phi_y.*repmat((TempD)',1,Q),1);

        % X_s
        X_sH = eqm.X_sH;

            % export
            B1 = sum(e^sigma*D_F./D1.*X1);  
                
        % D_s(q)
        D_s = phi_s.*X_sH/sum(phi_s.*Psig);

        % D(q) (new guess)
        D_Hnew = D_m+D_s;

        % c(q) (new guess)
        c_new = (newmbar.*sum(phi_y.*repmat(Psig,Q,1),2)').^(1/(1-sigma));   
        
        % update
        D_Hnew = D_H + 0.2*(D_Hnew - D_H);
        c_new = c + 0.2*(c_new - c);
        iter = iter+1;      
        
        % check convergence
        err1 = max(abs(D_Hnew-D_H)./D_H);
        err2 = max(abs(c_new-c)./c);
        err = err1 + err2;

    end
    
    % Warning
    if iter == maxIter
        msg = 'User Warning: Maximum Iterations Reached.';
        warning(msg)
    end 

    % Adjust back
    Pi0 = Pi0/rescale;
    Pi1 = Pi1/rescale; 
    
    % Save the Export Decision
    ExportCutoffQ = (k0*lnzQ-log(sigma*gamma)+log(kron(ones(Nf,1),Pi1)-kron(ones(Nf,1),Pi0))-mu_E)/sigma_E;
    ExportProbQ = normcdf(ExportCutoffQ);        
       

    
%% Return

    eqm_cf.Pi0 = Pi0;
    eqm_cf.Pi1 = Pi1;
    eqm_cf.D_H = D_H;
    eqm_cf.D_F = D_F;
    eqm_cf.D0 = D0;
    eqm_cf.D1 = D1;
    eqm_cf.D_s = D_s;
    eqm_cf.D_m = D_m;
    eqm_cf.c = c;
    eqm_cf.Psig = Psig;
    eqm_cf.P = Psig.^(1/(1-sigma));
    eqm_cf.X_m = X_m;
    eqm_cf.X0 = X0;
    eqm_cf.X1 = X1;
%     eqm_cf.V = V;
%     eqm_cf.Vbar = Vbar;
%     eqm_cf.M = M;
%     eqm_cf.xi = xi;
%     eqm_cf.theta_v = theta_v;
%     eqm_cf.theta_m = theta_m;  
%     eqm_cf.rv_ads = rv_ads;
%     eqm_cf.rv_rev = rv_ads.*D_H./D1; 
%     eqm_cf.C_v1 = C_v1;
    eqm_cf.B1 = B1;
    eqm_cf.P_s = P_s;
    eqm_cf.e = e;
    eqm_cf.X_s = X_sH;
    eqm_cf.X = sum(X0+X1); 
    eqm_cf.ExportCutoffQ = ExportCutoffQ;
    eqm_cf.ExportProbQ = ExportProbQ;
    eqm_cf.iter_inner = iter;
    eqm_cf.err_inner = err;
    
    eqm_cf.newmbar = newmbar;    
    
    F = eqm_cf;

    
end